library(ggplot2)
rm(list=ls())
sel=list.files(pattern=".snp")
sel1=sel
col=paste("total_reads>",5:15,sep="")
row=gsub(".snp","",as.character(sel1))
result=data.frame(matrix(NA,length(row),ncol=length(col)),row.names = row)
names(result)=col

for(i in 5:15){
file=read.table(sel1[1],head=T,sep="\t")
file=file[file$exits=="T",]
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
file$con_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)
file$case_reads=as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
#sel=which(file$con_reads>=i | file$case_reads>=i)
sel=which(file$total_reads>=i)
file=file[sel,]
result[1,i-4]=dim(file)[1]
vardata1=data.frame(normal_var_freq=file$normal_var_freq,type=c(rep("DC_unaffected",length(file$normal_var_freq))))
vardata2=data.frame(tumor_var_freq=file$tumor_var_freq,type=c(rep("DC_affected",length(file$normal_var_freq))))
p=ggplot(vardata2,aes(x = tumor_var_freq,colour=type))+
geom_line(stat = "density")+
geom_vline(xintercept = 0.25,lty=6)+
geom_vline(xintercept = 0.75,lty=6)
p=p+geom_line(data=vardata1,aes(x=normal_var_freq,colour=type),stat = "density")
for(j in 2:13){
file=read.table(sel1[j],head=T,sep="\t")
file=file[file$exits=="T",]
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
file$con_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)
file$case_reads=as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
#sel=which(file$con_reads>=i | file$case_reads>=i)
sel=which(file$total_reads>=i)
file=file[sel,]
result[j,i-4]=dim(file)[1]
vardata1=data.frame(normal_var_freq=file$normal_var_freq,type=c(rep("DC_unaffected",length(file$normal_var_freq))))
vardata2=data.frame(tumor_var_freq=file$tumor_var_freq,type=c(rep("DC_affected",length(file$normal_var_freq))))
p=p+geom_line(data=vardata1,aes(x=normal_var_freq,colour=type),stat = "density")
p=p+geom_line(data=vardata2,aes(x=tumor_var_freq,colour=type),stat = "density")
}
p=p+xlab("variant allele frequency")+theme_light()
fileName=paste("T_total_reads",i,".png",sep="")
png(file=fileName,width=700,height=650)
print(p)
dev.off()
}
write.csv(result,"data_log_T.csv",quote=F)


###
for(i in 5:15){
file=read.table(sel1[1],head=T,sep="\t")
file=file[file$exits=="T",]
file$tumor_var_freq=as.numeric(gsub("%","",file$tumor_var_freq))/100
file$normal_var_freq=as.numeric(gsub("%","",file$normal_var_freq))/100
file$total_reads=as.numeric(file$normal_reads1)+as.numeric(file$normal_reads2)+as.numeric(file$tumor_reads1)+as.numeric(file$tumor_reads2)
sel=which(file$total_reads>=i)
file=file[sel,]
result[1,i-4]=dim(file)[1]
}